1 



Heatbath Noise Methods in Lattice QCD 

Walter Wilcox a 

a Department of Physics, Baylor University, Waco, TX 76798-7316 



o 
o 

(N 
Or 



> 
oo 
O 

o 

On 
O 



X 
3 



In a recent paper, de Forcrand has pointed out that matrix inversions using Gaussian noise need not be iterated 
to full convergence, but instead may be solved approximately and treated as a heatbath. Gaussian noise however 
is not optimal for minimizing variance. It shown here how his algorithm may be generalized to a mixture of 
Gaussian and Z(N) noise, resulting in a smaller effective variance for some operators. 



1. Introduction 

Usually, one solves A(f> = ff, where ff is a noise 
vector, to full convergence. DeForcrand pointed 
outjl} that this is not necessary for Gaussian 
noise. Formally, 
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Introduce an auxilliary field 
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Consider a change, 

=t-$,A( = x-r,x = A$+ff, (3) 

where ff is a complex Gaussian noise vector and r 
is the residual vector in the solution for <j>' . One 
can show the change is accepted with probability 
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where 

AS = \A$'\ 2 + \x-A^\ 2 -\A^ 
= 2Re(r t ■ A{$ -$')). 
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With the assumption that A(f>, A(j>' , and fare 
uncorrelated Gaussian vectors of variance N, N 
and e 2 N respectively (N is the dimensionality of 
A), DeForcrand shows that 
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The computational overhead is sim- 



dot products per acceptance check. This can save 
a factor of 2 to 3 in computer time. This is the 
general idea; generalizations are also presented by 
DeForcrand. 

2. Accelerating Z(N) Noise 

Gaussian noise is not optimal for signal extrac- 
tion j2|. Therefore it is of interest to see if heat- 
bath methods can be adapted to use a mixture of 
Gaussian and Z(N) noise (Z(2) used here). 

For this purpose, we begin with the expression, 
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where Z % is a particular Z(2) noise vector and 
TV = 2 N is the number of Z(2) noises in the vector 
space and a 2 + b 2 = 1. Let O = J2 ^ A ^ 
can do the integrals to get, 
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Then with 
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we see that the answer is just TrA -1 , but with a 
weighting over Gaussian (a 2 ) and Z(2) (b 2 ) noises. 
We introduce x as before, 

< O >= , „L - / BfV / D$ 



ply one matrix- vector product plus several matrix 
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Figure 1. The acceptance, Pa, as a function of 
logio(e) for b/a = (boxes) and b/a = 1 (circles). 

treating Z % as a dynamical variable, the change 
in the action is now, 

AS = \[\A$ - bZ 1 ' | 2 + \x- A$' + bZ l '\ 2 
a z 

-\A$-bZ l \ 2 -\x- A$+bZ l \ 2 ]. (11) 
This is again a heatbath, with 

AS = 2Re(r« t • (A a $ - -Z l - A a $' + -Z 1 ')), (12) 
a a 

where A a = A/ a and r* = f/a. This is of the 
same form as above and has the same acceptance 
since (A a $-±Z l ) and (A a 4>' -^Z -1 ' ) have variance 
N. The rescaled r°" is used to define the residual, 
e, in the computer program. 

Although the acceptance is the same, the num- 
ber of iterations is greater for a given cutoff, e, on 
the new residual vector since it is defined by di- 
viding by a < 1. However, since the convergence 
on the residual is typically exponential, changes 
in a are accomodated by a modest number of ex- 
tra iterations. The fact that the acceptance is the 
same at the rescaled e value we will see has the 
helpful consequence that one does not have to re- 
search for the optimum e value at which to run, 
even though now we use a noise mixture. 

The above discussion can be generalized to jus- 
tify more complicated exponential shifts A<f> — > 
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Figure 2. The number of mixed Z(2)/Gaussian 
heatbath iterations divided by the model accep- 
tance, N mixed I Pa, as a function of e in units of 
10~ 5 for b/a = (boxes) and b/a = 1 (circles). 



A<f) — C % where C % can depend on other noises, 
parameters, etc., besides the specific one C l = 
bZ l used above. I am using a solver which has 
enforced even/odd preconditioning in the Wil- 
son/Dirac matrix M. One can show that in this 
case one is not simulating A = M directly as in 
the above discussion, but a different, shifted sys- 
tem and the residual vector is purely in one sector 
(even/odd) or the other. 

Figs. 1 and 2 show acceptance data on a 12 3 x24 
lattice using the Wilson/Dirac matrix at k = 
0.148. I am calculating the average acceptance 
on 100 noises at various cutoff values of e and 
with b/a = and 1. For this solver I get about 
a 6% increase in the number of iterations for this 
ratio. (The model suggests about 14% increase 
in the number of iterations for an 80% mixture 
(b/a = 2) and about a 30% increase for a 97.3% 
mixture (b/a = 6) for the parameters in this sim- 
ulation.) Fig. 1 shows Pa in the crucial window 
of e of from 10 -4 through 10~ 3 , where it varies 
between ~ 0.97 to 0.65 for b/a — 0,1 for a single 
gauge configuration. The acceptance is the same 
for both b/a values within errors and agrees with 
Eq. (g) . The reason this interval is crucial is illus- 
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trated in Fig. 2. There we see that the number of 
mixed Z(2)/Gaussian heatbath iterations divided 
by the acceptance, N m i xed / Pa, has a minimum 
at e ~ 4 x 10~ 4 for b/a = 0, 1. Shown here also is 
the model value for this quantity given by 

N mlxed /P A = iV , (13) 

erfc(eyf ) 

where N m i xed is assumed determined by \\f\\ = 
r e^ cNm "" i . As noted above, I am using a solver 
which has enforced even/odd preconditioning in 
the Wilson/Dirac matrix M. The upshot for this 
simulation is that TV in ([|) and (^) must include 
a factor of |, N = | X (12 3 x 24 x 4 x 3). 

Fig. 3 shows the normed effective mixed 
Z(2)/Gaussian variance in the operator ipiji as a 
function of b/a. It has a very shallow minimum at 
b/a w 10. By normed effective variance I mean 
the ratio {V mixed x N mixed ) / (V Gauss x N Gauss ), 
which takes into account that N m i xed > Ncauss 
for b 0. The model and data suggest that the 
minimum of this ratio is « 0.3. The horizontal 
line gives the fully converged value of the variance 
ratio Vz(2) /Vcauss, which is sa 0.2. Thus, one 
looses a factor of « 2/3 in the effective variance 
compared to the fully converged Z(2) simulation. 
The gain in computer time is reduced from a fac- 
tor of from 2 to 3 to a factor of from 1.33 to 2. 
These numbers are apparently typical. Other op- 
erators with larger converged Vz{i)/Vgo,uss ratios 
have sharper minima at smaller b/a. 

3. Conclusions and Observations 

Ref. jlj shows that heatbath methods can speed 
up simulations of many disconnected loop opera- 
tors or by a factor of « 2 to 3. However, Gaussian 
noise is not optimal and so heatbath methods do 
not help for operators whose variance is diago- 
nally dominate, such as Wilson ipip. It has been 
shown here that these methods can be generalized 
to a mixture of Gaussian and Z(N) noise. With 
an mixture/iteration penalty factor of about 2/3, 
the Gaussian noise accelerates the Z(N) sector. 

For diagonally dominate operators, there exists 
a Z(N)/Gaussian ratio that minimizes the vari- 
ance. The noise can then be tailored to the opera- 




Figure 3. Effective variance of the 6^0 simula- 
tion for ipip normalized to the Gaussian simula- 
tion with b = as a function of b/a. 

tor ("designer noise"). The optimum b/a param- 
eter can be numerically estimated from the pa- 
rameters of the model and an independent mea- 
surement of the fully converged Z(N) /Gaussian 
variance ratio. 

It is not shown here, but the even/odd struc- 
ture of the Wilson/Dirac matrix can be exploited 
to increase the computer time gain of diago- 
nally dominate operators further by restricting 
the Gaussian noise to only one sector or the other. 
This and other aspects of heatbath noise methods 
will be discussed in a future publication. 
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